// File:       complex.c++
// Version:    1.00
// Author:     (c) Miles Sabin, 1997
// Purpose:    approximation to ANSI C++ complex template

// Some of the following is due to AT&T and hence carries the
// following copyright notice and disclaimer.

// Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
//
// Permission to use, copy, modify, and distribute this software
// and its documentation for any purpose and without fee is hereby
// granted, provided that the above copyright notice appear in all
// copies and that both that the copyright notice and this
// permission notice and warranty disclaimer appear in supporting
// documentation, and that the names of AT&T Bell Laboratories or
// Bellcore or any of their entities not be used in advertising or
// publicity pertaining to distribution of the software without
// specific, written prior permission.
//
// AT&T and Bellcore disclaim all warranties with regard to this
// software, including all implied warranties of merchantability
// and fitness.  In no event shall AT&T or Bellcore be liable for
// any special, indirect or consequential damages or any damages
// whatsoever resulting from loss of use, data or profits, whether
// in an action of contract, negligence or other tortious action,
// arising out of or in connection with the use or performance of
// this software.

// Change log:
//  06/04/97   v. 1.00

#include "complex.h"

#include "exception.h"
#include "iostream.h"
#include "sstream.h"


// Implementation of complex<T>

template<class T>
complex<T>& complex<T>::operator*=(complex<float> const& rhs)
  {
    T temp = (re_*rhs.real())-(im_*rhs.imag());
    im_ = (re_*rhs.imag())+(im_*rhs.real());
    re_ = temp;
    return *this;
  }

template<class T>
complex<T>& complex<T>::operator*=(complex<double> const& rhs)
  {
    T temp = (re_*rhs.real())-(im_*rhs.imag());
    im_ = (re_*rhs.imag())+(im_*rhs.real());
    re_ = temp;
    return *this;
  }

template<class T>
complex<T>& complex<T>::operator/=(complex<float> const& rhs)
  {
    T ar = fabs(rhs.real());
    T ai = fabs(rhs.imag());
    T nr;
    T ni;
    T t;
    T d;

    if(ar <= ai)
    {
      t = rhs.real()/rhs.imag();
      d = rhs.imag()*(1+t*t);
      nr = (re_*t+im_)/d;
      ni = (im_*t-re_)/d;
    }
    else
    {
      t = rhs.imag()/rhs.real();
      d = rhs.real()*(1+t*t);
      nr = (re_+im_*t)/d;
      ni = (im_-re_*t)/d;
    }

    re_ = nr;
    im_ = ni;

    return *this;
  }

template<class T>
complex<T>& complex<T>::operator/=(complex<double> const& rhs)
  {
    T ar = fabs(rhs.real());
    T ai = fabs(rhs.imag());
    T nr;
    T ni;
    T t;
    T d;

    if(ar <= ai)
    {
      t = rhs.real()/rhs.imag();
      d = rhs.imag()*(1+t*t);
      nr = (re_*t+im_)/d;
      ni = (im_*t-re_)/d;
    }
    else
    {
      t = rhs.imag()/rhs.real();
      d = rhs.real()*(1+t*t);
      nr = (re_+im_*t)/d;
      ni = (im_-re_*t)/d;
    }

    re_ = nr;
    im_ = ni;

    return *this;
  }


// Implementation of complex<T> free fns

template<class T>
complex<T> operator/(T const& x, complex<T> const& y)
{
  T ar = fabs(y.real());
  T ai = fabs(y.imag());
  T nr;
  T ni;
  T t;
  T d;

  if(ar <= ai)
  {
    t = y.real()/y.imag();
    d = y.imag()*(1+t*t);
    nr = x*t/d;
    ni = -x/d;
  }
  else
  {
    t = y.imag()/y.real();
    d = y.real()*(1+t*t);
    nr = x/d;
    ni = -x*t/d;
  }

  return complex<T>(nr, ni);
}

template<class T>
basic_istream_char& operator>>(basic_istream_char& is, complex<T>& x)
{
  basic_istream_char_sentry sentry(is);
  DESTROY_ON_THROW(sentry);

  if(sentry)
  {
    char c = 0;

    if(is.peek () == '(')
      is >> c;

    T re;
    T im = 0;

    is >> re;

    if(c == '(')
    {
      is >> c;
      if(c == ',')
        is >> im >> c;
    }

    if(c != 0 && c != ')')
      is.setstate(ios_base::failbit);
    else if(is.good())
      x = complex<T>(re, im);
  }

  return is;
}

template<class T>
basic_ostream_char& operator<<(basic_ostream_char& os, complex<T> const& x)
{
  basic_ostringstream_char s;
  s.flags(os.flags());
  s.precision(os.precision());
  s << '(' << x.real() << "," << x.imag() << ')' << ends;
  return os << s.str();
}

template<class T>
complex<T> pow(complex<T> const& x, int y)
{
  if(y == 0)
    return complex<T>(1.0);

  complex<T> r(1);
  complex<T> x_(x);

  if(y < 0)
  {
    y = -y;
    x_ = T(1)/x_;
  }

  while(true)
  {
    if((y&1) != 0)
      r *= x_;

    if((y >>= 1) != 0)
      x_ *= x_;
    else
      return r;
  }
}

template<class T>
complex<T> sqrt(complex<T> const& x)
{
  T r = abs(x);
  T nr;
  T ni;

  if(r == 0)
  {
    nr = r;
    ni = r;
  }
  else if(x.real() > 0)
  {
    nr = sqrt(0.5*(r+x.real()));
    ni = x.imag()/nr/2;
  }
  else
  {
    ni = sqrt(0.5*(r-x.real()));

    if(x.imag() < 0)
      ni = -ni;

    nr = x.imag()/ni/2;
  }

  return complex<T>(nr, ni);
}
